Elevating open data: 3D visuals of NYC open climate data with elevatr & rayshader

nyc_border <- st_read("https://data.cityofnewyork.us/resource/wh2p-dxnf.geojson")
inundation_border <- st_read("https://data.cityofnewyork.us/resource/5xsi-dfpx.geojson")

# add topography limited to borders of city and inundation zone
create_topo <- function(st, bbox, z_value) {
  elev <- get_elev_raster(st, z = z_value) 
  elev_masked <- mask(elev, bbox) # remove NAs to reduce matrix dims
  st_matrix <- raster_to_matrix(elev_masked)
  return(st_matrix)
}

nyc_matrix <- create_topo(nyc_border, nyc_border, 10)
inundation_matrix <- create_topo(inundation_border, inundation_border, 10)

zscale <- 10
n_frames <- 150
inundation_waterdepths <-  max(inundation_matrix, na.rm = TRUE)/2 - max(inundation_matrix, na.rm = TRUE)/2 * cos(seq(0,2*pi,length.out = n_frames))
img_frames <- paste0("surge", seq_len(n_frames), ".png")
for (i in seq_len(n_frames)) {
  message(paste(" - image", i, "of", n_frames))
  nyc_matrix |>
    sphere_shade(texture = create_texture("#b7f4f5", "#4e3e42","#8a9299","#195f60","#c2d1ce"), zscale = zscale) |>
    add_water(detect_water(nyc_matrix, cutoff = 1.3)) |>
    add_shadow(ray_shade(nyc_matrix, zscale = zscale)) |>
    plot_3d(heightmap = nyc_matrix,
            zscale = zscale,
            solid = FALSE,
            shadow = TRUE,
            water = TRUE,
            watercolor =  "#3289a0",
            theta = -5,
            phi = 65,
            zoom = 0.4,
            windowsize = c(1200, 900)
            )
   render_water(inundation_matrix, waterdepth = inundation_waterdepths[i]/zscale,
               watercolor ="#3289a0", zscale = zscale, remove_water = FALSE)
  render_snapshot(img_frames[i])
  rgl::clear3d()
}

magick::image_write_gif(magick::image_read(img_frames), 
                       path = "nyc_inundation.gif", 
                       delay = 5/n_frames)

nyc_border <- st_read("https://data.cityofnewyork.us/resource/wh2p-dxnf.geojson")
inundation_border <- st_read("https://data.cityofnewyork.us/resource/5xsi-dfpx.geojson")

# create bounding box to focus on the beach
riis_coords <- c(-73.893936,40.557853,-73.857586,40.578100)
names(riis_coords) = c("xmin", "ymin", "xmax", "ymax")
riis_bbox <- st_bbox(riis_coords)
riis_polygon <- st_as_sf(st_as_sfc(riis_bbox), crs = 4326)

# add topography limited to borders of beach and inundation zone
create_topo <- function(st, bbox, z_value) {
  elev <- get_elev_raster(st,  z = z_value, src = "gl1") # requires OpenTopography API key (https://opentopography.org/developers)
  elev <- crop(elev, bbox)
  elev_masked <- mask(elev, bbox) # remove NAs to reduce matrix dims
  st_matrix <- raster_to_matrix(elev_masked)
  return(st_matrix)
}

# higher zoom value = higher resolution
riis_matrix <- create_topo(nyc_border, riis_polygon, 14)
inundation_matrix <- create_topo(inundation_border, riis_polygon, 14)


zscale <- 10
n_frames <- 75
inundation_waterdepths <-  max(inundation_matrix, na.rm = TRUE) - max(inundation_matrix, na.rm = TRUE) * cos(seq(0,2*pi,length.out = n_frames))
img_frames <- paste0("surge", seq_len(n_frames), ".png")

for (i in seq_len(n_frames)) {
  message(paste(" - image", i, "of", n_frames))
  
  riis_matrix |> 
    sphere_shade(texture = "imhof3", zscale = zscale) |> 
    add_water(detect_water(riis_matrix, cutoff = 1.3)) |>
    add_shadow(ray_shade(riis_matrix, zscale = zscale)) |> 
    plot_3d(heightmap = riis_matrix, 
            zscale = zscale, 
            solid = FALSE,
            shadow = FALSE,
            water = TRUE,
            watercolor =  "#3289a0",
            theta = -5, 
            phi = 65
    )
render_water(inundation_matrix, waterdepth = inundation_waterdepths[i]/zscale,  
    watercolor ="#3289a0", zscale = zscale, remove_water = FALSE)
render_snapshot(img_frames[i])
rgl::clear3d()
}

magick::image_write_gif(magick::image_read(img_frames), 
                       path = "riis_inundation.gif", 
                       delay = 5/n_frames)

file.remove(list.files(pattern = "surge"))